unit Arenstorf;
interface
uses Matrix, ODEModel, Math;

type TArenstorf = class(TODEModel)
public
  function F(t: extended; y:MatrixPtr): MatrixPtr; override;
  constructor Create(init:MatrixPtr; problemSize:integer);
end;

implementation

  constructor TArenstorf.Create(init:MatrixPtr; problemSize:integer);
  begin
    Inherited Create(init, problemSize);
  end;

  function TArenstorf.F(t: extended; y:MatrixPtr): MatrixPtr;
  var
  values : MatrixPtr;
  mu, mu_, D1, D2 : extended;
  begin

    mu  := 0.012277471;
    mu_ := 1.0 - mu;
    D1  := Power( sqr(GetMatrixElement(y, 3, 1)  + mu)   + sqr(GetMatrixElement(y, 4, 1)), 3.0/2.0);
    D2  := Power( sqr(GetMatrixElement(y, 3, 1)  - mu_)  + sqr(GetMatrixElement(y, 4, 1)), 3.0/2.0);

    values := CreateMatrix(4, 1);

    SetMatrixElement(values, 1, 1, GetMatrixElement(y, 3, 1) + 2.0*GetMatrixElement(y, 2, 1) - mu_*(GetMatrixElement(y, 3, 1)+mu)/D1 - mu*(GetMatrixElement(y, 3, 1)-mu_)/D2);
    SetMatrixElement(values, 2, 1, GetMatrixElement(y, 4, 1) - 2.0*GetMatrixElement(y, 01, 1) - mu_*GetMatrixElement(y, 4, 1)/D1 - mu*GetMatrixElement(y, 4, 1)/D2);
    SetMatrixElement(values, 3, 1, GetMatrixElement(y, 1, 1));
    SetMatrixElement(values, 4, 1, GetMatrixElement(y, 2, 1));
    result := values;
  end;
end.
 